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An idea for evaluating transition probabilities in chemical reaction systems is proposed, which is 
efficient for repeated calculations with various rate constants. The idea is based on duality relations; 
instead of direct time-evolutions of the original reaction system, the dual process is dealt with. 
Usually, if one changes rate constants of the original reaction system, the direct time-evolutions 
should be performed again, using the new rate constants. On the other hands, only one solution 
of an extended dual process can be re-used to calculate the transition probabilities for various rate 
constant cases. The idea is demonstrated in a parameter estimation problem for the Lotka-Volterra 
system. 


I. INTRODUCTION 

Recent developments of experimental techniques en¬ 
able us to obtain detailed data for bio-chemical reactions 
within cells, and the importance of the role of data anal¬ 
ysis has been increased. In small systems such as cells, it 
has been shown that discrete characteristics could play 
important roles [l| , and hence chemical reactions in such 
small systems should be treated as chemical master equa¬ 
tions. That is, conventional rate equations are not ade¬ 
quate, in which the noise effects are neglected and num¬ 
bers of chemical substances are approximated as contin¬ 
uous variables. Hence, it would be necessary to treat 
discrete variables directly. The time-evolution of the dis¬ 
crete variables is directly expressed via the chemical mas¬ 
ter equations. In the past, various analytical and numeri¬ 
cal methods for the chemical master equations have been 
developed (for example, see 0.) 

Transition probabilities in the chemical master equa¬ 
tions are one of the important quantities, especially in 
time-series data analysis. For example, consider an esti¬ 
mation problem for rate constants from time-series data 
of chemical reactions. When all time-series data for num¬ 
bers of chemical substances and chemical reactions are 
available, the Bayesian statistics gives easily the estima¬ 
tion of the rate constants However, it is in general dif¬ 
ficult to obtain such detailed time-series data, and partial 
observations with discrete time steps should be consid¬ 
ered in realistic experiments. In general, the partial ob¬ 
servation cases need many numerical time-evolutions for 
different rate constants, and large computational costs 
are needed. For example, in [^, an estimation method 
based on gradient descent techniques has been proposed, 
and the method needs many iterated calculations to seek 
rate constants with the largest likelihood value. Al¬ 
though sometimes the estimation problems have been 
performed via approximations [a-lTj, these approxima¬ 
tions treat discrete variables as continuous ones, and then 
the approximations could be inadequate for small sys¬ 
tems such as cells. Hence, it is useful to develop more 
rapid and direct evaluation methods for transition prob¬ 
abilities for different rate constant cases. 

In the present paper, an idea to evaluate transition 


probabilities in chemical master equations is proposed. 
The idea is based on duality concepts; instead of direct 
time-evolutions for the original chemical master equa¬ 
tions, the corresponding dual processes are evaluated. 
Employing extensions of states, it is possible to construct 
an extended dual process which does not explicitly in¬ 
clude some rate constants of the original processes. Fur¬ 
thermore, it will be shown that only a time-evolution for 
the extended dual process enables us to evaluate transi¬ 
tion probabilities for the original chemical master equa¬ 
tions with various rate constants. This characteristics 
of the extended dual processes could reduce largely the 
computational costs for the calculations of the transition 
probabilities for various rate constants. 

This paper is organized as follows. In Sec. [Ill the basic 
idea for the usage of the dual process is explained. Sec¬ 
tion imi gives a demonstration for the derivation of the 
dual processes using the Lotka-Volterra system. Con¬ 
cluding remarks are given in Sec. IIVI 


II. DUALITY RELATIONS FOR TRANSITION 
PROBABILITIES 

In this section, the idea based on the duality relations 
is explained. A concrete example for the derivations of 
the dual process will be given in Sec. IHII 


A. Stochastic chemical reaction systems 

Consider a stochastic chemical reaction system with 
M chemical substances, Xi,... ,Xm- Denote n = 
(ni,..., tim) as the state of the reaction systems, where 
rim denotes the number of chemical substance Xm- 
Note that rim takes any positive integer or zero {rim € 
{0,1,2,•••}.) Assume that there are R chemical reac- 


2 


tions, which take the form 

uiiXi + • ■ • uimXm —^ + • • • vimXm, 

M21-^1 + ■ ■ • U2mXm > V 21 X 1 H V2mXm, 


know transition probabilities for various rate constant 
cases. Hence, in the following discussions, the chemical 
master equations are investigated from the viewpoint of 
bosonic operators, which enables us to obtain an idea to 
avoid the repeated time-integrations. 


URiXi + ■ ■ ■ UrmXm —^ VRiXi + ■ ■ ■ VrmXm, 

( 1 ) 


B. Doi-Peliti formalism 


where Urm and Vrm correspond to the stoichiometries as¬ 
sociated with the m-th reactant and product of the r-th 
reaction, respectively. The r-th chemical reaction has a 
rate constant Cr- The actual rate for the r-th reaction de¬ 
pends on the number of chemical substances, and a rate 
law or hazard, hr{n,Cr), is introduced Q. For example, 
for 2Xj -I- Xk A Xi, the hazard should be defined as 


h{n,c) = 



and so on. In the present pa¬ 


per, the following redefined hazard, h'^{n), is introduced: 


hr{n, Cr) = Crh'^(n), 


( 2 ) 


where 




= Cr 


M 

n 


(3) 


and 


K(n) 


M 


n 


rim'- 


(4) 


These notations are convenient to describe the extension 
of the dual process later. In addition, using the quantities 
{urrn} and {vrm}^ the net effect reaction matrix A is 
defined Q; the components of H, {arm}, are defined as 

^rm — '^rm '^rm- 

Using the above notations, the chemical master equa¬ 
tions are written as follows: 




R 

= ^ [Crh'r{n — Ar)P{n 
r—1 


Ar,t) 


Crh'r{n)P{nA)] , 

(5) 


where Ar is the r-th row of the matrix A. For details of 
the master equations, see, for example, ii- 

When one employs the direct numerical time- 
integration for the original chemical master equations 
in Eq. @ (with a suitable truncation for finite num¬ 
bers of equations), the transition probability for a cer¬ 
tain initial state n to final state m with fixed rate con¬ 
stants {ci,...,c/j} can be evaluated. If one wants to 
know the transition probabilities for different rate con¬ 
stant cases, the numerical time-integrations must be per¬ 
formed again with the new rate constants; these repeated 
time-integrations are time-consuming when one wants to 


The chemical master equations in Eq. m are essen¬ 
tially the infinite number of coupled ordinary differential 
equations. There are some analytical methods to rewrite 
the chemical master equations, and one of the methods is 
the so-called Doi-Peliti formalism [8l-[l^. The Doi-Peliti 
formalism has been widely used, ranging from the re¬ 
search fields of reaction-diffusion systems m to genetic 
switches [3 [ 3 . The method is based on the algebraic 
probability theory 03, and the following bosonic cre¬ 
ation operators for the m-th chemical substance, and 
annihilation operators, am, are used: 


a-ln] 


= 1, 

(6) 

aln] 

— [^m 5 — O5 


(7) 

0-ln' 

] ~ [^m, Rm'] 



■ [am, am'] — [Rm,R\n'] 

= 0, for m ^ mb 

(8) 


That is, the creation and annihilation operators for the 
same chemical substance Xm do not commute with each 
other; for different chemical substances, these operators 
can commute. 

The actions of the creation and annihilation operators 
on state |n) = |ni,..., um) in a Fock space are defined 
as 

1 ^1, ... , Rm, ■ • ■ , Rm) — \Ri , ■ ■ ■ , Rm T 1, ■ . . , RM ), 


(9) 

am \ R1, ' • ' , Rm, ■ • ■ , RM ) — Rm 1^1 , ... , Rm , Rm) ■ 

(10) 

The corresponding dual (bra) states {m\ = (mi,..., trmI 
are introduced as satisfying the following inner product: 

{m\n) = 6{ni,Rii) ■ ■ ■ 6{nM,RT-M)n'., (11) 

where S(rim,Rim) is the Kronecker delta, and 

n! = m! ■ • • Rm'-- (12) 


It is easy to confirm that the actions of the creation and 
annihilation operators to the bra states become as fol¬ 
lows: 


(ni,.. 

■ 5 5 • • 

• I'^M l^m — • • • t 15 ■ ■ - ; | ? 



(13) 

(ni,.. 

■ 5 5 • • 

■ ; l^m — (^1 j • • ■ ; 15 ■ ■ • : \ ■ 



(14) 
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Using the above notations, the state \P{t)) is defined 
as 

oo oo 

l-PW) = XI ■■■ X! (15) 

m—O riM—O 

In order to derive the time-evolntion equation for the 
state \P{t)), the following quantities are introduced, 
which correspond to the redefined hazards in the orig¬ 
inal chemical master equations: 

M 

nr=l[ , (16) 

m—1 

M 

n (17) 

m—1 


In addition, noting (rici | = (0|(®ci)"°U we have 

{nci\zi) = (24) 

Second, the following time-evolution operator U®’' is 
introduced instead of the original C: 


^ex _ 




(25) 


where 

CT = «ci (Hi - . (26) 

Furthermore, the state \P{t)) is extended as 

OO OO 

\P^^{t))= X ••• X PM\n)\z,), (27) 

ni =0 riM—O 


Then, the original chemical master equations in Eq. m 
are rewritten as follows: 

||P(t)) = £|P(t)), (18) 


C = J2^r, (19) 

r— 1 

and 

Cr = Cr {Pr - . ( 20 ) 

Finally, the transition probability Pn^m{t) is written 
in terms of the Doi-Peliti formalism as 

Pn^mit) = -^{m\P{t)) = —{m\e‘^*\n), (21) 

m! m! 

where e^*\n) corresponds to a solution of the time- 
evolution starting from state n. Note that we need the 
factor 1/m! because of the characteristics of the inner 
product in Eq. m- 


C. Extended dual process 

Here, additional bosonic operators are introduced to 
obtain an extended dual process, in which some rate con¬ 
stants are vanished and additional states emerge. For 
simplicity, here only one rate constant ci will be replaced 
with the additional bosonic operator; extensions for mul¬ 
tiple cases are straightforward. 

First, additional bosonic operators, and are 
introduced. Here, the following property of the coherent 
states is important: 

ac^\zi) = zi\zi), (22) 

where \zi) is the coherent state with parameter zi G R, 
which is defined as 


(23) 


and the time-evolution for the extended state \P™{t)) 
obeys 

^|P®’'(t))=r’^|P'=’'(t)). (28) 

Note that U®’' in Eq. (E51) gives the same quantity with C 
in Eq. (fT8l) when zi = ci, because of the characteristics 
of the coherent state in Eq. (1221) . 

Third, the following bra state is defined: 

OO OO OO 

Wt)\= X X'" X <(>(«’ »^ci,i)(nci I (n|, (29) 

nc;^=0ni=0 riM—O 


where (n| = (ni,... ,nM\, and (n-cj corresponds to the 
state related to the additional bosonic creation and an¬ 
nihilation operators and 0 ^. 

Fourth, instead of the time-evolution for |P®’'(t)), the 
following time-evolution for the bra state is considered: 

|(</(^)| = (</(t)|£®^ (30) 

that is, 

||^(t)) = (rX|</(t)), (31) 

where (P®’')^ is the conjugate of U®’'. Equation (1^ is 
written in terms of the Doi-Peliti formalism, and it is 
possible to derive the corresponding infinite coupled or¬ 
dinary differential equations for {^(n, ricj, t)}. That is, 
introducing the following ‘hazard’: 


M 


Kin) = 


i'^m '^rm)- 


(32) 


we have 
—0(n,nci,t) 

= h[(n + Ar)(/>(n -|- Ar^Uc^ — l,t) — /i(.(n)(/(n, rici “ lU) 
R 

+ X [crKin + Ar)4>in + t ^Ci ; t) - CrK{n)4){n,nc-,,t) 

r^2 

(33) 
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Here, note that the time-evolution equations for 
rici, t)} 9-'^® n.ot the chemical master equations in 
general; the equations do not satisfy the probability con¬ 
servation law, and then {</>(«, rici, t)} cannot be inter¬ 
preted as a probability distribution. Although it is pos¬ 
sible to recover the probabilistic nature by using similar 
discussions used in [l^, it is enough to use 71 ^, t)} 

to calculate the transition probabilities here. 

The main idea in the present paper is the replacement 
of the time-evolution for \P{t)) with {(j>{t)\. Hence, 

= = 0|(m|e^'**‘|n)|zi = ci) 

m! 

^ 00 00 00 

= XI X'" X </'(«',i)Ki|(n'|n)|zi = Cl) 

ric^—On'^—O n'j^—0 
I 00 

= ^ X (34) 

Kci =0 

where the initial condition for the bra state is 

, . . . , , TT-cj , t = 0) 

_|l for n'l = mi,... ,n'^^ = mM,nci = 0, 

10 otherwise. 

The above discussions imply the following fact: The 
transition probabilities with the rate constant ci can be 
evaluated from Eq. (IMll using the solutions of the ex¬ 
tended dual process. Note that the solutions for the ex¬ 
tended dual process, {</>(«, rici, t)}, does not include the 
rate constant ci explicitly. Hence, only one numerical 
time-integration for the extended dual process is neces¬ 
sary to evaluate the transition probabilities n —> m for 
various ci cases. 


III. DEMONSTRATION OF THE DERIVATION 
AND APPLICATIONS OF DUAL PROCESSES 


Here, a demonstration for the derivation of the ex¬ 
tended dual process is shown by using the famous Lotka- 
Volterra system, which has been already used for param¬ 
eter estimation problem in Q: 


'Xi^2Xi, 

A 1 +A 2 4 ^ 2 , 
Xi+X2-^Xi + 2 X 2 , 
X2^%. 


(36) 


The chemical master equation for the Lotka-Volterra 


system is written as 

dP{ni,n2,t) 

dt 

= a [{ni - l)P[ni - 1, n 2 , t) - niP(ni, 712 , t)] 

-I- P [(tIi -I- 1 ) 712^(711 -I- l,7l2,t) - 7ll7l2H(7ll,7l2,t)] 

+ 5 [ 711(712 — l)P(7ll,7l2 — l,t) — 771712^(711,712, t)] 

+ l[{n2 + 1)P(77i, 772 -b l,t) - 772^(771,772, t)] , (37) 

where 77i and 772 are the numbers of particles Xi and 
A 2 , respectively. The time-evolution operator in the Doi- 
Peliti formalism is defined as 


C =a{a\a\ai — a\ai) -I- f3{aia\a2 — a|aiaja 2 ) 

+ d(a|aia2a2a2 — a|aia2a2) + 7(02 — a|a2). ( 38 ) 

We assume that precise values for the rate constants a 
and P are unknown. Therefore, the two bosonic operators 
(oa and ap) and their adjoint operators (aj, and a^) are 
introduced. Using certain constants a' and /?', a and P 
in Eq. (|38)) are replaced as a'aa and P'ap, respectively, 
and then the following extended time-evolution operator 
is derived: 

=a'aa{a\a\ai — ajoi) -|- P'ap[aia\a 2 — a\aia\a 2 ) 
-b i5(a|aia2a2a2 — alaiala2) -b 7(02 — a\a2)- 

( 39 ) 

Note that the replacement a a'ua is used here, instead 
of a —>■ Oq, (and the same as P). These replacements 
do not mean approximations; they correspond to simple 
variable transformations. The reasons to introduce the 
replacements are as follows: 

• Sometimes we know the rough values (or only the 
orders) of unknown rate constants; these additional 
information can be embedded into the extended 
dual process with small modifications of the dis¬ 
cussions in Sec. H, as we will see here. 

• As we will see below (Eq. (l44l) l. the final expression 
corresponds to the Taylor-type (Maclaurin-type) 
expansion. Hence, in practical, it is important to 
use small aja' and PjP' to confirm the rapid con¬ 
vergence of the summation in Eq. (I44|l . Hence, 
sometimes the replacements reduce the practical 
computational issues. 


The adjoint operator for is 



and using the discussions in Sec. mi the extended dual 
process obeys the following time-evolution equation: 


||</)(t)) = (£-)n<^(t)), 


(41) 
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FIG. 1. Measurement data. Symbols denote partial obser¬ 
vations with discrete time steps Ar°^'* = 0.3. For compari¬ 
son, the original processes are plotted as dotted lines. This 
artificial time-series data is generated with the parameters 
a = 0.2, P = 0.02, 5 = 0.02, and 7 = 0.3. 


where 

\m 

00 00 CXD CXD 

= E E E E (I>[ni,n2,na,ni3,t)\ni)\n2)\na)\ni3). 

na =0 n/ 3=0 ni =0 n 2—0 

(42) 

Then, the following coupled ordinary differential equa¬ 
tions are derived: 

d 

= a'[(ni -I- l)ni(/)(ni -|- l,n 2 ,nQ - 
- ni(j){ni,n 2 , - 1, np, t)] 

+ - 1 , n2, na,np- 1 , t) 

- nin24>{ni,n2,na,np - l,t)] 

+ S[ni{n2 + I)n24>{ni,n2 + 1 , na,np,t) 

- nin2(/)(ni,n2,na,n/3,t)] 

-I- 7 [(j){ni,n 2 - l,na,np,t) - n 2 (j){ni,n 2 ,na,np,t)]. 

(43) 

Finally, the solutions of Eq. (gSl) are used to evaluate a 
transition probability from state ni,n 2 to state mi, m 2 
as follows: 


P, 


ni ,712—7712 

ni!n2! 


it) 


00 00 / \ 

i: i: (7)”" 7 

1 —n -rj n \ ^ / 


, y 

mi!TO2! 

nQ=0 7ljp—0 


710 


0 (^ 1 7 ^2 7 7 ^/3 7 t ) . 

(44) 


Notice that Eq. (H41) has a form of the Taylor-type 
expansion around the origin. Hence, the usage of the 


dual process corresponds to calculations of the expan¬ 
sion coefficients. In addition, it is easy to calculate the 
derivatives with respect to the parameters from the same 
(l>{ni,n 2 ,na,np, t); there is no need to perform additional 
time-evolution. 

The derived formula for the transition probabilities in 
Eq. (|44ll can be, for example, used in parameter esti¬ 
mation problems. Assume that we have discrete-time 
observations of time-series date, which are depicted in 
Fig. [TJ The observation-time interval is = 0.3, 

and the full time-series data are not available. That is, 
only the observation values and at discrete time 
I = 1, 2,..., N} are available, where = 

At°^^ = 0.3, and TV = 31 in Fig. [T] 

Assume that two parameters S and 7 are known. The 
other two parameters, a and /3, should be estimated from 
the time-series data. In order to seek the parameters, the 
following likelihood function is calculated; 

N-l 

^c(,p) = Prob(n^*+^yy^y^*\ 4 *y,/?), ( 45 ) 

i=l 

where Prob(nJ*^^\ a, (?) is the probabil¬ 

ity of the state change with parameters a and /?. Ad¬ 
ditionally, it is usual to calculate the following log- 
likelihood function, 

i{a,P) = log C{a,P), (46) 

instead of the original likelihood function. 

In general, the log-likelihood function in Eq. (H51) 
should be re-calculated for various values of a and /3, and 
the tasks need high computational costs. In contrast, the 
formula in Eq. (IT41) is suitable to depict the contour or 
heatmap for the log-likelihood function. That is, values 
of the log-likelihood function can be evaluated only from 
a time-integration of the extended dual process. Hence, 
there is no need to repeat the time-integration for differ¬ 
ent parameters. (Here, assume that the orders of scales of 
a and /3 are previously known, and the settings a' = 1.0 
and P' = 0.1 are used in Eq. (l43)l .'l Additionally, the 
information of the derivatives of the log-likelihood func¬ 
tions is easily obtained from Eq. (|44)) ; there is no need to 
perform additional time-integrations. 

In order to demonstrate the usage of the duality rela¬ 
tions, the numerical calculations for heatmaps and null- 
clines for the log-likelihood functions are performed as 
follows. Both Eqs. (1571) and (H51) are coupled ordinary dif¬ 
ferential equations, and the numbers of the equations are, 
in principle, infinite. Usually, the time-integration needs 
a finite cut-off for states; here, only states with 0 ^ 69 
are considered for each particle. I checked that the finite 
cut-off is enough for the truncation of the summation in 
Eq. (H41) . For the time-integration for the dual processes, 
the usual 4th-order Runge-Kutta method is employed. 
Figures Ha) and (b) are the results for the heatmap; 
Fig. Ht>) is the enlarged one of Fig. Ha)- Figure H^) 
shows the directions of derivatives of the log-likelihood 
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(b) 



- 122.02 

-122.08 

-122.14 

- 122.20 

-122.26 

-122.32 

-122.38 

-122.44 

-122.50 


(c) 0.026 

0.025 

0.024 

0.023 

0.022 

0.24 0.25 0.26 0.27 0.28 0.: 

a 



a 

t t I > > t t i 1 i i t • . . ■ • -^- r-*-^ 

1 t t t t I I » » ! 1 ^ - : 1 yK i i 

i t t t t J » i i . . ; . . t t i * * 

i t t I I 1 1 1 ; ; . : t I y * I * * 

i t t I : 1 ^ ; . i i t t y i i i i I 

\ \ I \ ^ \ ^ y t * I / ^ 

■ t t j y i i I \ 

- ^ i i i i i y i I /^ ^ ^ 
i i i i i t iyf t I \ \ \ \ \ 

! i i i I I I \X i \ \ \ ^ t t 

. 1 I t ; i I I \y I ^ t t < < ♦ < t 


■ t * i I ♦ 1 i I '’A^ ^ t t < 1 ^ 1 

• t t i 1 » * 1 i 

i t I i ; ; i ttitiiit-- 

i t i V V t t t ♦ i t t ; 

i i i t t t ♦ I t t : - • ; • 1 

\ f i/T i t t t t t •••<111 

>y t t \A \ \ 1 t \ ■ ■■■■■■■ \ t \ 

t i i A i \ \ ‘ - ';;.*• ‘ t » i i t 

t i A t t • ' - . t t f » i < 

f t ^ ‘^-i ! i t t t t > i 1 


FIG. 2. Heatmaps and nullclines for the log-likelihood func¬ 
tion. (a) and (b) show the heatmaps. (b) is an enlarged one; 
note that the color profiles in (a) and (b) are different, (c) 
shows the directions of derivatives and nullclines. 


FIG. 3. Plots of the quantities defined in Eqs. (EU and 
(1521) (up to n = 30.) When numerical evaluation for larger 
n cases is performed, the intercept with 1/n = 0 gives the 
radius of convergence. The plots correspond to the cases with 
m = 37, n 2 = 11, mi = 36, m 2 = 12. For bI^\ j3' = 0.1 and 
P = 0.0236 are used; a' = 1.0 and a — 0.261 are used for 
rW 

J->n 


coefficients: 


(a) ni\n 2 \ ^ 


mi\m2\ 


miimo! Vq / 


UQ—O 


Ud—O 


(47) 

(48) 


Note that and depend on ni, n 2 , m-i, m 2 , /3, P' 
and t; the dependencies are not explicitly shown in 
and for notational simplicity. Then, the transition 
probability is written in the following power-series with 
one variable: 


P, 


n\ ,n2—^m\ ,m2 


II 

/ a \ ^ 

Vav 

(49) 

00 

n—O 

y- 

(50) 


In order to discuss the radius of convergence of the 
power-series, one may use the so-called Domb-Sykes plot 
[ 1 ^ . Here, we use the method introduced by Mercer 
and Roberts M because complicated patterns appear in 
{cn^'^} and In the Mercer-Roberts method, the 

following quantities are calculated: 


functions, and nullclines. Again, note that there is no 
need to perform additional time-integrations here; only 
the same solutions for the extended dual process are re¬ 
peatedly re-used to depict all Figs.[2Ka), (b), and (c). 

Finally, I give some comments for the radius of con¬ 
vergence of the power series in Eq. (ITi)) . Equation (H41) 
has two variables, and we here introduce the following 






(a) (a) 

C „_2 ■ 

J/3) J/3) 

^n+l^n — 1 


SP)SP) 

^ ^n-2 


- {Xr 
AX 

- {ciX 

iPP) f 2 ’ 


(51) 

(52) 


for n = 2 ,3,4,.... As discussed in [T^, for example, the 
reciprocal of the radius of convergence, l/r, for Eq. 
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is given as the intercept with 1/n = 0 when we plot 
versus 1/n. Figure [3] shows the analysis; the plots corre¬ 
spond to the case with rii = 37, n 2 = 11, mi = 36, m 2 = 
12, a' = 1.0, a = 0.261, /3' = 0.1, and /3 = 0.0236. 
(Other parameters give the similar behaviors.) The re¬ 
sults in Fig. [3] would not be enough to obtain the precise 
values of the radius of convergence; the coefficients with 
larger n in the power-series are needed, but there are 
limitations of the memory capacity in the Runge-Kutta 
methods and numerical precision. However, as shown in 
Fig. [31 the intercepts with 1/n = 0 seem to take near¬ 
zero values, and the radius of convergence would be large 
enough. 


IV. CONCLUDING REMARKS 

As shown in the present paper, the duality relations en¬ 
able us to reduce the repeated time-integrations for var¬ 
ious rate constants. This feature is applicable, for exam¬ 
ple, to obtain heatmaps for the log-likelihood functions; 
compared to the direct time-integrations for the original 
system, only one time-integration for the extended dual 
process is enough, as demonstrated in the present paper. 
In addition, the derivatives of the log-likelihood functions 
can also be obtained easily by using the same numerical 
solutions for the extended dual process. Of course, it is 
also possible to calculate the Jacobian matrix, the Hes¬ 
sian matrix, and so on. 

It should be noted here that the number of random 
variables in the extended dual processes is larger than the 
original one. Hence, the computational cost for one time- 
integration becomes larger than the original one. How¬ 
ever, if we can perform the time-integration for the ex¬ 
tended dual process with reasonable computational costs, 
the numerical results can be repeatedly re-used, which 
will finally reduce the whole computational costs. 

Here, the current limitations of the usage of the duality 
relations should be stated. 

1. As discussed above, the usage of the duality re¬ 
lations corresponds to the Taylor-type expansion. 


Hence, we must pay attention to the convergence. 
In practice, it is preferable that the expanded vari¬ 
ables take values smaller than 1. In order to avoid 
this convergence issue, it is useful to use the re¬ 
placements (variable transformations) in Sec. HI. 

2. In the demonstration, the coupled ordinary dif¬ 
ferential equations for the dual process should be 
solved numerically. When the numbers of param¬ 
eters in the expansion are large, it becomes im¬ 
possible to solve the coupled ordinary differential 
equations (the curse of dimensionality.) Hence, at 
this stage, the current approach is suitable to see 
the behavior of the log-likelihood functions with a 
few varying variables. 

3. In order to treat many variable cases, one may won¬ 

der if the Monte Carlo approach is available. It 
is, in principle, true; although the time-evolution 
equation for the dual process does not correspond 
to a stochastic process in general, the stochas¬ 
tic nature can be recovered by using more exten¬ 
sions (see [ 13 .) Actually, as for the duality be¬ 
tween the stochastic differential equations and the 
birth-death processes, the Monte Carlo approach 
has already been employed However, the 

Taylor-type expansion needs sometimes solutions 
with high precisions, and the Monte Carlo approach 
is still time-consuming. The usage of the impor¬ 
tance sampling may avoid this problem; this is be¬ 
yond of the scope of the present paper, and under 
investigation. 

At the present moment, the duality relations are use¬ 
ful to investigate cases with a few varying variables. In 
future works, it is important to develop approximation 
methods or efficient numerical methods for the extended 
dual processes. 
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